library(glmnet)
library(survival)
library(tidyr)
library(dplyr)
library(ggplot2)
inputdata<-read.csv("select_features.csv",header = T, check.names = F)


x<-as.matrix(inputdata[,c(4:55)])



y<-data.matrix(Surv(inputdata$OS,inputdata$Event))
#lasso回归
lasso<- glmnet(x, y, 
               family = 'cox', 
               nlambda=1000, 
               alpha=1)
#10折交叉验证
lassoCV <- cv.glmnet(x, y, family = "cox",
                     type.measure = "deviance",
                     nfolds =10)


par(mfrow=c(1,2))
plot(lasso, xvar="lambda")
plot(lassoCV)

se_ambda<-lassoCV$lambda.min

se_coef<-coef(lassoCV,s="lambda.min")

se_coef
index<-which(se_coef!=0)

coef<-se_coef[index]

diffvariables = row.names ( se_coef )[index] # 
lasso.result.se <- cbind(diffvariables , coef)
lasso.result.se<- as.data.frame ( lasso.result.se )
lasso.result.se$coef <- as.numeric(lasso.result.se$coef)


ggplot ( aes ( x = reorder ( diffvariables , coef ), y = coef , fill = diffvariables ), data = lasso.result.se )+
  geom_col ()+
  coord_flip ()+
  theme_bw ()+
  labs ( x ="")+
  ggtitle (" LASSO identified variables ")+
  scale_fill_brewer ( palette ="Set1")+
  theme ( legend.position = "")





# 提取plot(fit.lasso)中的数据
#提取lasso中的lambda和回归系数并建立数据框
plot.data <- data.frame(as.matrix(lasso$lambda), as.matrix(t(lasso$beta)))
#将款数据格式转换为长数据格式以便作图
plot.data%>%
  rename(lambda='as.matrix.lasso.lambda.')%>%
  pivot_longer(cols = 2:ncol(plot.data),
               names_to ="variable",
               values_to = "coef")->plot.data
#这样我们就得到了绘图数据，有三列，分别为lambda、每个变量和他们的回归系数



# 绘制图形
library(viridis)
# 使用viridis调色板
colors <- turbo(57)#我们有76个变量
p4<-ggplot(plot.data, aes(x = log(lambda), y = coef, color = variable)) +
  geom_line(size=1,alpha=0.8) +
  scale_color_manual(values =colors)+
  #scale_x_log10(label=scientific_10)+
  labs(title = "LASSO Regression Path",
       x = "Log lambda",
       y = "Coefficient") +
  theme_bw()+
  theme(legend.position = "")
p4



# 提取cv.glmnet数据
cv.df <- data.frame(lambda =lassoCV$lambda,#交叉验证中的lambda
                    mse =lassoCV$cvm,#mse
                    sd =lassoCV$cvsd)#sd


# 绘制10折交叉验证过程图
p5<-ggplot(cv.df, aes(log(lambda), mse)) +
  geom_point() +#点图
  geom_errorbar(aes(ymin = mse - sd, ymax = mse + sd), width = 0.1) +#添加误差棒
  scale_x_continuous(name = "Log lambda") +
  scale_y_continuous(name = "Mean Squared Error") +
  ggtitle("10-fold Cross-validation using Lasso Regression")+
  geom_vline(xintercept = log(lassoCV$lambda.min), linetype = "dashed", color = "#E41A1C",size=1) +
  geom_vline(xintercept = log(lassoCV$lambda.1se), linetype = "dashed", color = "#377EB8",size=1) +
  annotate(geom = "text",label=("lambda.min: 16 non-zero variables"),x=-7.5,y=12,color = "#E41A1C",size=3)+
  annotate(geom = "text",label=("lambda.lse: 1 non-zero variables"),x=-7.5,y=11.9,color = "#377EB8",size=3)+
  theme_bw()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))
p5 



library(ggrisk)
library(survival)
library(rms)

#colnames(LIRI)

LIRI<-inputdata[,c("OS","Event", "DL-1" , "DL-5",  "DL-54" ,"DL-78" ,"DL-82")]
colnames(LIRI)<-gsub("-","_",colnames(LIRI))

fit <- cph(Surv(OS,Event)~DL_1+DL_5+DL_54+DL_78+DL_82, LIRI)
ggrisk(fit,cutoff.value='cutoff',
       color.A=c(low='#00BFFF',high='#D2691E'),#A图中点的颜色
       color.B=c(code.0='#00BFFF',code.1='#D2691E'), #B图中点的颜色
       color.C=c(low="#00BFFF",cutoff='#FFFAF0',high='#D2691E'))

surv.dat<-LIRI

f = survival::coxph(Surv(OS, Event) ~ ., data = surv.dat)
#  2.risk point nomgram.points and lp
riskscore = f$linear.predictors




##TCGA

library(survival)
library(survivalROC)
library(survminer)



riskscore <- data.frame(riskscore = as.numeric(riskscore),
                        group = ifelse(riskscore >1.2,"HRisk","LRisk"), # 根据中位数分组
                        row.names = rownames(surv.dat),
                        OS.time = surv.dat$OS,
                        OS = surv.dat$Event,
                        stringsAsFactors = F)

Sur <- Surv(riskscore[,3] ,riskscore[,4])
sfit <- survfit(Sur ~ riskscore[,2],data=as.data.frame( riskscore[,2]) )
ggsurvplot(sfit,
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T,
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_risk','Low_risk'))+
  labs(x = "Days")

###HX_cohort


LIRI<-read.csv("HX_DLfeature.csv",header = T, check.names = F)
colnames(LIRI)<-gsub("-","_",colnames(LIRI))

fit <- cph(Surv(OS,Event)~DL_1+DL_5+DL_54+DL_78+DL_82, LIRI)
ggrisk(fit,cutoff.value='median',
       color.A=c(low='#00BFFF',high='#D2691E'),#A图中点的颜色
       color.B=c(code.0='#00BFFF',code.1='#D2691E'), #B图中点的颜色
       color.C=c(low="#00BFFF",median='#FFFAF0',high='#D2691E'))

surv.dat<-LIRI

f = survival::coxph(Surv(OS, Event) ~ ., data = surv.dat)
#  2.risk point nomgram.points and lp
riskscore = f$linear.predictors






library(survival)
library(survivalROC)
library(survminer)



riskscore <- data.frame(riskscore = as.numeric(riskscore),
                        #group = ifelse(riskscore >1.2,"HRisk","LRisk"), # 根据中位数分组
                        row.names = rownames(surv.dat),
                        OS.time = surv.dat$OS,
                        OS = surv.dat$Event,
                        stringsAsFactors = F)

res.cut <- surv_cutpoint(riskscore, time = "OS.time",
                         event = "OS",
                         variables = names(riskscore)[1],
                         minprop = 0.3)

res.cat <- surv_categorize(res.cut)
riskscore$group<-res.cat$riskscore

Sur <- Surv(riskscore[,2]*30 ,riskscore[,3])
sfit <- survfit(Sur ~ riskscore[,4],data=as.data.frame( riskscore[,4]) )
ggsurvplot(sfit,
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T,
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_risk','Low_risk'))+
  labs(x = "Days")
